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ABSTRACT: Stress and strain fields in a two-dimensional pixelwise disordered 
system are computed by a Fast Fourier Transform method. The system, a model for a 
ductile damaged medium, consists of an elastic-perfectly matrix containing void pix- 
els. Its behavior is investigated under equibiaxial or shear loading. We monitor the 
evolution with loading of plastically deformed zones, and we exhibit a nucleation / 
growth / coalescence scenario of the latter. Identification of plastic "clusters" is eased 
by using a discrete Green function implementing equilibrium and continuity at the 
level of one pixel. Observed morphological regimes are put into correspondence with 
some features of the macroscopic stress / strain curves. 
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1 INTRODUCTION 

The versatile Fast Fourier Transform (FFT) method of Moulinec, Suquet, and Michel 
represented a breakthrough in the computation of the stress and strain fields in linear or 
nonlinear composites. This method uses the Lippmann-Sch winger integral equation 
of the strain field in a homogeneous reference linear medium, written in the Fourier 
space. Nonlinearity is embedded via a pointwise heterogeneous polarization field 
which depends on the constitutive law, and which is computed in direct space. The 
integral equation is solved iteratively through FFT and inverse FFT transformations 
for which efficient routines are available. The ill-convergence of the basic iterative 
procedure is alleviated by use of an "augmented Lagrangian method" and by Uzawa's 
algorithm, see [1] for details. 

With some modifications, this method is used hereafter to investigate the build- 
up and incipient localization of plastic deformation in elastic-perfectly plastic porous 
pixelwise disordered systems [2]. Such systems consist of an on-lattice realization of 
a random system where the material properties of adjacent material elements are sta- 
tistically uncorrelated from point to point, in the limit where the size of the material 
element goes to zero. Other realizations of this type of disorder (in the bond form) 
include random spring or resistor networks. Such systems are particularly attractive 
as benchmarks for homogenization methods, since their microcopic disorder correla- 
tion length is the smallest possible. Then, the self-consistent linear effective medium 
approximation is exact to fourth order in correlations (in the sense of diagram ex- 
pansions, e.g. 0]), as has been shown for dielectric media or for random resistors 
networks (RRNs) Q. 



2 DISCRETE GREEN FUNCTION IN FFT CALCULATIONS 



A two-dimensional (2D) square lattice of pixels (i, j) of size L 2 is considered. Pixels 
are either randomly chosen as voids, in concentration /, or as elastic-perfectly plastic 



matter elements with flow stress Y. Deformation theory [4] is used. Using the con- 
tinuum Green tensor of the strain (GT) as in fll] (we call this approach "CGI"), we 
observe that : (i) convergence of the FFT method is slow for infinite contrast; (ii) the 
fields, and the plastically deformed zone where the Mises norm of the stress cr eq is 
locally equal to the flow stress Y exhibit a spurious "checkerboard" pattern (see Fig. 
[U which complicates a subsequent identification of plastic clusters (cf. Sec. 3). 

Discrete Green functions — In the alternative solution considered here, introducing 
the unit vectors along the axes e 1 such that = Sij (the Kronecker symbol), we 
follow the RRN scheme of Ref . Q] and use forward and backward finite-difference 
schemes for the compatibility and equilibrium equations (discrete Fourier Transforms 
in elasticity are also used in Ref. 0]): 

^•(x) = - [^-(x + e*) - Uj{x) + (i ^ j)] , t^iW ~ a ^( x ~ &)] = 0. 

3 

(i) 

Contrary to a centered-difference scheme, equilibrium is enforced here on any finite 
connex subset of pixels. In the discrete Fourier representation with discrete Fourier 
momenta qi = 2-KmijL, ra* = 0, 1, . . . L — 1 (the conventions of Ref. [3] are used), 
these equations read [here i = (— l) 1 / 2 ]: 

e k i = (i/2) (k k ui + kiu k ) , ikfa k i = 0, with kj = 2 sm(qj/2)e iqj /2 , (2) 
where * denotes the complex conjugate. The associated GT for the strain is: 

G ijkl (q) = - (N^hK + N^kjk* k + N7 k %kt + N-'kjkt) /4 (3) 

where the acoustic tensor is Nij = k% Cikijh with Cijki the elastic tensor of the 
isotropic reference medium of Lame moduli fi and A. The restriction to generalized 
plane strain loading (fields independent of z) is obtained by setting q% = and retain- 
ing only the indices i,j,k,l = 1,2. Though an analytic expression of G is available 
(2[]> the inversion of N is most easily carried out numerically. The continuum GT 
(e.g. (J]) is retrieved in the long wavelength limit #2 <^ 1- A centered-difference 
scheme would instead lead to a real GT, simply obtainable from the continuum GT by 
replacing its Fourier momenta by kj = sin(^j). 

The GT ([3]) does not comply with the square symmetry of the grid, and lacks 
major symmetry, due to its nonzero imaginary part: G^ki — G\ u - . The inconsistency 
comes from the fact that the equilibrium equation in (Q]) is a balance condition for 
forces transmitted by bonds linking nearest-neighbor pixels, whereas the constitutive 
law used is appropriate for material points only, i.e. the pixels. As a result, the forward 
and backward directions on a each cartesian axis are not equivalent, which leads to 
asymmetric field patterns. Two easy work-arounds are considered. The first one, 
called discrete Green 1 (DG1), consists in replacing © by its symmetrized version 

G i]li {Gijki(qi,q2) + Gijki(-qi, -#2) 

+ [G ijkl (- qi ,q 2 ) + G ijkl (qi,-q 2 )] } /4. (4) 

This GT corresponds to no discretization scheme in direct space, but should be inter- 
preted as the GT of some non-local medium in which the stiffness tensor C^ki is re- 
placed by a non-local convolution kernel with finite range of order 1 . The second way, 
called DG2, consists in carrying out four different calculations on the same system up 
to final convergence, employing each one of the four Gijki (#1, #2), Gijki (—91, —^2), 
Gijki {— Qi ? #2), Gijki (Qi 7 —#2) in turn, an in taking the average of the four converged 




Figure 1: Differences between the methods. Periodic case with one void under 
equibiaxial loading. The gray field shows the equivalent strain field norm, e eq = 

[(2/3)e^ -e^ -] 1 / 2 , where e'- = eij — 5ijSkk/3. Pixels near the void surface with the 

highest field values have been thresholded out (~ 1.5% of the total number). The right 
image is an enlargement. 



strain fields as the final result. We also consider, for comparison purposes, field pic- 
tures obtained from CGI by means of a local average post-processing consisting in 
taking 5 -point averages on the current pixel and its four nearest neighbors. We call 
this the CG2 method. 

Setting k 2 = k\ + k 2 , the 2D displacement field is obtained from £$j(q) using 
(the mode q = is irrelevant): 

% % 
^i(q) = (*i [ £ n ~ ^22] + 2/c 2 £i2} , u 2 (q) = {k 2 [en - £22] - 2ki e 12 } • 

In the computations hereafter, the bulk and shear moduli are K = 1 and /1 = 0.4 and 
Y = 0.5. The loading is prescribed by imposing an overall strain, which is increased 
until the stress reaches its flow value [Q]]. 

Differences between the methods are illustrated in Fig. [T] where the equivalent 
shear strain in a periodic medium under equi-biaxial loading is displayed, in quadrants 
of a unit cell with one circular void (volume fraction / = 0.1) Q Method CG2 does not 
suppress the "checkerboard" pattern of method CGI. Method DGl blurs excessively 
the strain field. The best result without "checkerboard" effect is obtained with DG2. 
The strain field inside the void depends on the method, but is physically irrelevant. 

Convergence issues — Two convergence indicators are used. First, we require the 
stress divergence (computed in Fourier representation) to be such that (| | diver 1 1 2 ) < 
r]f (a) : (a). Additional steps of the iterative algorithm are then carried out until 
(<r n+1 - a n ) : (tr n+1 - a n ) < j]\ (a n+1 ) : (tr n+1 ) between steps n and n + 1. 
The overall prescribed tolerance is specified by the pair (771, 772). For benchmarking, 
typical values of 771,2 of order 10 -5 , and system sizes L = 512, 1024, 2048 were 
considered. For both the periodic void lattice and the pixelwise disordered medium 
on which we focus hereafter, method DG2 converges faster (i.e., in fewer iterations) 
than CGI. The absolute precision on the average stress, extrapolated to infinite system 
sizes by an inverse power law fit of the size dependence, is then typically 10 -3 . For 

lr The strain localization motif, of a surprisingly rich sub-structure - remark the bands of finite width 
which bound the localization zone - is approximately made of logarithmic spirals leaving the void surface 
at angle 45° |4] as predicted by slip-line theory (e.g., U), though the lines are somehow distorted and 
blurred by lattice effects. 



more precise calculations, values of 771,2 of order 10 8 and system size L = 4096 
have been used (see GJ for details). 



3 BUILD-UP OF PLASTIC DEFORMATION IN THE POROUS MEDIUM 

Using method DG2, we investigate the link between the overall stress-strain curve and 
the development of the plastic zones where a eq = Y in a pixelwise disordered porous 
medium. The medium can be considered as made of three phases: elastic (where 
cr eq < Y), porous (where a = 0), and plastic (cr eq = Y). Under increasing loading, 
plastic zones develop starting from the voids (around which the shear stress is larger), 
grow, coalesce and eventually "percolate" in the system. Void or plastic clusters, 
are identified (with the Hoshen-Kopelman algorithm [6]) as pixel sets of same phase 
connected by the nearest-neighbour criterion. The "void-plastic" phase comprises the 
voids and the plastic zones. In this phase the local tangent shear modulus is zero. 
Due to local unloading effects (possibly an artefact of the - reversible - deformation 
theory), purely plastic clusters in the vicinity of the pores can temporarily disconnect 
during the first loading stages, especially under shear loading. This has bearing on 
cluster counting. This was corrected by modifying the first-neighbor connectivity rule: 
we furthermore connect to the nearest void, each plastic cluster not already connected 
to a void. Each cluster then contains one void at least. 
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Figure 2: Equi-biaxial loading. Left: stress-strain curve and geometric indicators 
as a function of the applied overall strain. Right: empirical stress-loading formula 
compared to FFT stress/strain curve. 

Equi-biaxial loading — Fig. 2 (left) displays, for a porosity / = 0.01, various geo- 
metrical indicators, as a function of the applied equi-biaxial overall strain sq = (e m ), 
along with the macroscopic equibiaxial stress/strain curve cjq = (cr m )(eo) (the brack- 
ets denote a spatial average), and the opposite of its second derivative —d 2 G^jde\\ 
the normalized volume fraction of void-plastic zone f v = v p /(l — /) where v p is 
the volume fraction of void-plastic zone; the volume fraction v max of the largest void- 
plastic cluster; the proportion of isolated voids clusters i(eo) (i.e. voids not connected 
to a plastic zone); the number n(so) of void-plastic clusters (normalized to the num- 
ber of voids at so = 0); and the coalescence rate r(so) = dn(eo)/deo, multiplied 
by a magnifying factor so as to make it conspicuous on the curve. Four characteristic 
regimes are isolated, marked on the figure: (1) elastic loading regime: the number 
of isolated voids remains constant; (2) growth of plastic zones around the voids: the 
number of isolated voids diminishes, / p grows approximately quadratically, the coa- 
lescence rate develops a huge peak (coalescence between plastic clusters originating 
from neighboring voids). The regime ends up as i(eo) = 0, all voids having devel- 
oped a plastic zone; (3) regime of stabilized coalescence: f p increases linearly, the 
coalescence rate somewhat stabilizes, then decreases more slowly as a largest void- 



Figure 3: Equi-biaxial loading for porosity / = 0.1. Left: largest void-plastic cluster 
(white) at percolation (so — 0.9). Right: equivalent norm of the shear strain field at a 
larger overall strain (so = 1.1). 



plastic cluster emerges; (4) stress saturation regime: it begins at the percolation of the 
void-plastic zone where the largest plastic cluster grows fast, and where / p increases 
slower than linearly with the strain. From these observations, we see that the second 
derivative of the stress-strain curve is strongly correlated to the coalescence rate, and 
that the macroscopic flow stress attains its final order of magnitude at percolation of 
the plastic zone through the medium. Typical plastic clusters and shear strain field are 
displayed in Fig. [3] Their fractal character will be discussed elsewhere. The shear 
strain is strongly localized (but not uniform) within the plastic zones. 




Figure 4: Simple shear loading. Same legend as fig. 2. 



Loading in simple shear — Fig. 4 (left) displays the same quantities as Fig. 2, for the 
same medium, in simple shear loading where £o — (e xy ) and ao — {o'xyY The main 
difference is the disappearance of the stabilized coalescence regime. Fig. \5\ displays a 
typical instance of the fields for a weak porosity, with straight shear bands. 

Empirical formula for the stress-strain curve — Amusingly, one can reproduce the 
stress-strain curves using the above geometric indicators and the porosity-dependent 

effective compressibility and shear elastic moduli K and /L We indeed arrived at 
the empirical formulae (the effective elastic moduli are computed on the simulated 
system): 

a ~ K(f) £ ck { + « [« P (eo) - v™(e )] } , (5) 



Figure 5: Simple shear loading at / = 10 -4 . System size L = 1024 at applied strain 
£o = 2. Left: plastic (white) and elastic (black) zones. Right: equivalent norm of the 
shear strain field. 



where © and © apply to pressure and shear loadings respectively, and where a is a 
fitting number different in each case. In both expressions, the term containing cluster 
numbers reproduces the zone of maximal curvature of the stress-strain curve (regime 
of maximal coalescence rate), whereas the plastic volume fraction enters the descrip- 
tion of the saturation regime. These formulae are compared to the FFT stress / strain 
curves in Figs. 2 and 3 (right). Though the present "guesswork" should not be taken 
too seriously (in particular, the dependence of a with respect to / has not been studied, 
and the plastic volume fraction does not enter both formulas in the same manner), such 
an approach might nonetheless ultimately contribute to enrich the effective-medium 
approach to porous media. Indeed taking correlations into account so as to incorporate 
micro structural information in the effective-medium framework is a notoriously diffi- 
cult problem. The possibility of using global indicators of geometric nature such as 
the fraction of plastic zone, or cluster numbers, for which phenomenological evolution 
models could be proposed might then constitute a useful alternative. 

Acknowledgments: Y.-P.P. thanks Pierre Suquet for drawing his attention to Ref. O. 
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